package scu.maqiang.possion;

import scu.maqiang.csparsej.tdouble.Dcs_cholsol;
import scu.maqiang.csparsej.tdouble.Dcs_common;
import scu.maqiang.csparsej.tdouble.Dcs_compress;
import scu.maqiang.fes.*;
import scu.maqiang.numeric.*;
import scu.maqiang.mesh.*;

public class CoefPossion2DT3 {
	public static void main(String[] args) {
		long startTime = System.nanoTime();
		Mesh2T3 Th = new Mesh2T3().square2D(200, 200, T3Type.Left);
		//
		// //Th.rotate(Math.PI/6);
		// //AbstractVectorFunction ff = new ComplexFunction();
		// //Th.scale(1, Math.PI * 3 / 4);
		// //Th.transform(ff);
		// //Th.writeToTecplot("mesh.dat");
		//
		FES2T31 fs = new FES2T31(Th);
		System.out.println("Number Of DOF " + fs.GetNdof());

		ScalarFunc rf = (xy, label, param) -> 5 + Math.sin(4 * Math.PI * xy[0]) + Math.sin(4 * Math.PI * xy[1]);
		double[][] coef = new double[1][];
		coef[0] = fs.valueFromFunc(rf, null);
		SRMatrix A = new SRMatrix(fs.GetNdof());
		fs.assembleStiff(coef, BVPType.COMMON, A);
		// Console.WriteLine("Number of A " + A.getNumber());
		//
		// ConstFunction source = new ConstFunction(10);
		double[] f1 = new double[fs.GetNdof()];
		fs.assembleSource(new double[]{10.0}, BVPType.COMMON, f1);
		// System.out.println("Size Of f1 = " + MVector.toString(f1));
		//
		// Func g = new FuncConst(1.0);
		// double[] flux = fespace.assembleHeatFlux(g, 1, QuadType.Gauss1);
		//
		double[] RHS = f1; // MVector.add(f1, flux);
		fs.applyBCZero(A, RHS, Direct.X, BCType.MBN, 1, 2, 3, 4);

		// fespace.applyDirichletHeat(A, RHS, f, 2);
		// fespace.applyDirichletHeat(A, RHS, f, 3);
		// fespace.applyDirichletHeat(A, RHS, f, 4);

		// System.out.println("A = " + MMatrix.toString(A.full()));


//		A.sort();
//		System.out.println(A.display());
//		NewIterSSolver solver = new NewIterSSolver(A, 1e-6, 5000);
//		double[] x0 = new double[fs.GetNdof()];
//		// solver.precondition(A, RHS);
//		solver.PCGSSOR(RHS, x0, 1.5, 1);
//		System.out.println(Arrays.toString(x0));
//
//		Th.toTecplot("scPossion.dat", x0);
//
//		System.out.println("min: " + MVO.min(x0) + "\tmax: " + MVO.max(x0));
		// fespace.writeToTecplot("solution.dat", x1);
		// System.out.println("The Norm = " + MVector.L2Norm(x0));
		// System.out.println("The Error = " + MVector.distance(x0, x1));
		// System.out.println("The RError is " + MVector.distance(x0, x1) /
		// MVector.L2Norm(x0));

		//A.sort();
		//System.out.println(A.display());
		Dcs_common.Dcs myDcs = A.toDcs_matrix();
		//Dcs_print.cs_print(myDcs, false);
		//Dcs_lusol.cs_lusol(0, myDcs, RHS, 1.0);
		Dcs_common.Dcs T = Dcs_compress.cs_compress(myDcs); /* A = compressed-column form of T */
		//System.out.print("T:\n");
		//Dcs_print.cs_print(T, false); /* print A */
		Dcs_cholsol.cs_cholsol(1, T, RHS);
		//System.out.println(Arrays.toString(RHS));
		System.out.println("min: " + MVO.min(RHS) + "\tmax: " + MVO.max(RHS));
		long endTime = System.nanoTime();
		System.out.println("Elapse Time: " + (endTime - startTime) / 1e9 + " s");
	}
}
